Recognizing communities among microbes, meiofauna and macrofauna at hydrothermal vents
Sheryl Murdock
25 November 2020
1 Study description
Microbial and faunal assemblages associated with the tubeworm Ridgeia piscesae were collected from thirteen tubeworm bushes in the Endeavour Hydrothermal Vent Marine Protected Area and Middle Valley aboard the CCGS Tully using the remotely-operated vehicle ROPOS in August of 2016. Macro- and meiofaunal species were identified and enumerated from 9 assemblage (“grab”) samples. Microbial community DNA was extracted from biofilm/detrital matter associated with all 13 assemblages. Microbial biomass was separated into two size fractions (20-64µm, 0.2-20µm) by serial filtration and two method replicates for each size fraction comparing biomass removal methods were included from assemblages EMw1 and EMw6. Microbial DNA was extracted from 7 diffuse fluid and 3 background fluid samples for a total of 40 DNA extractions.
2 Data Import
2.1 Faunal species data
The faunal dataset contains counts of individuals for 41 macrofaunal taxa (1mm sieve) and 25 mieofaunal taxa (64µm sieve; inlcudes juveniles of 11 macrofauna). Large samples were split, counted and scaled up to full sample numbers. Details of counted splits are in Supplementary Tables 1 and 2. Data are provided as counts for meiofauna and relative abundance for macrofauna due to planned inclusion of macrofaunal species counts in another manuscript. Please contact R. Boschen-Rose (reboschen-rose@outlook.com) for macrofaunal species counts.
# Macrofauna species counts (available by request)
macro9 <- read.csv("Macro9counts.csv", header = TRUE, row.names = 1)
macro9 <- as.data.frame(t(macro9))
# Meiofauna species counts
meio9 <- read.csv("Meio9counts.csv", header = TRUE, row.names = 1)
meio9 <- as.data.frame(t(meio9))2.2 Microbial OTU data
The microbial dataset contains OTU counts from reads produced on Illumina MiSeq runs using 16S/18S rRNA gene primers for Bacteria, Archaea, and microeukaryotes. Size fractions and method replicates are combined here. See “SamplingBiasTests” for justification.
# Bacterial OTU counts
codaBact23nonsingleton <- read.csv("BacterialOTUs.csv", header = TRUE,
row.names = 1)
codaBact23nonsingleton <- as.data.frame(t(codaBact23nonsingleton))
# Eukaryal OTU counts
codaEuk23nonsingleton <- read.csv("EukaryalOTUs.csv", header = TRUE, row.names = 1)
codaEuk23nonsingleton <- as.data.frame(t(codaEuk23nonsingleton))
# Archaeal OTU counts (ECw11 not included- too few reads)
codaArch22nonsingleton <- read.csv("ArchaealOTUs.csv", header = TRUE, row.names = 1)
codaArch22nonsingleton <- as.data.frame(t(codaArch22nonsingleton))Files with size fractions and replicates kept separate (40 individual DNA extracts). Singleton OTUs have been removed.
# Bacteria
Bact40 <- read.csv("Bact40nonsingleton.csv", header = TRUE, row.names = 1)
# Eukarya
Euk40 <- read.csv("Euk40nonsingleton.csv", header = TRUE, row.names = 1)
# Archaea (two extracts of ECw11 removed- too few reads)
Arch38 <- read.csv("Arch38nonsingleton.csv", header = TRUE, row.names = 1)QPCR-balanced microbial community. Three microbial domains were balanced to whole community relative abundances based on quantitative PCR results and multiplied by 105 to scale to whole number pseudo counts.
balMic <- read.csv("QPCR_balanced_microbes.csv", header = TRUE, row.names = 1)
balMic <- as.data.frame(t(balMic))2.3 Metadata
Sample_type codes: f=fluids, B=background, bk=cell removal through agitation in bucket of artificial seawater, d=direct removal of cells in biofilms. T_category for diffuse fluids is determined by the basal temperature of the associated assemblage, NOT the diffuse fluid temperature itself.
metadata23 <- read.csv("metadata23.csv", header = TRUE, row.names = 1)
metadata23## Source Sample_type Vent_field Temperature T_category Substratum
## EMw1 Tubeworms bk Main_Endeavour 33.2 high sulfide
## EMw2 Tubeworms d Main_Endeavour 5.3 low sulfide
## EMw3 Tubeworms bk Main_Endeavour 3.4 low basalt
## EMw4 Tubeworms d Main_Endeavour 37.1 high sulfide
## EMw5 Tubeworms bk Main_Endeavour 13.4 low sulfide
## EMw6 Tubeworms bk Main_Endeavour 2.8 low basalt
## EMw7 Tubeworms d Main_Endeavour 33.0 high sulfide
## EMw8 Tubeworms d Main_Endeavour 37.4 high sulfide
## ECw9 Tubeworms bk Clam_Bed 15.0 low sulfide
## ECw10 Tubeworms bk Clam_Bed 31.0 high sulfide
## ECw11 Tubeworms d Clam_Bed 5.2 low basalt
## MVw12 Tubeworms bk Middle_Valley 27.4 high sulfide
## MVw13 Tubeworms bk Middle_Valley 21.0 low sulfide
## EMf1 Diffuse f Main_Endeavour 5.0 low sulfide
## EMf2 Diffuse f Main_Endeavour 4.6 low sulfide
## EMf3 Diffuse f Main_Endeavour 4.1 low basalt
## EMf7 Diffuse f Main_Endeavour 12.0 low sulfide
## ECf9 Diffuse f Clam_Bed 2.6 low sulfide
## ECf11 Diffuse f Clam_Bed 6.1 low basalt
## MVf12 Diffuse f Middle_Valley 23.0 low sulfide
## bkEC Background B Clam_Bed 2.0 low none
## bkEM Background B Main_Endeavour 2.0 low none
## bkMV Background B Middle_Valley 2.0 low none2.4 Load necessary R packages
library(vegan)
library(reshape2)
library(ggplot2)
library(ggrepel)
library(ggpubr)
library(funrar)
library(compositions)
library(zCompositions)
library(propr)
library(ALDEx2)3 Sample diversity
Diversity was calculated using the Inverse Simpson metric. For microbial domains from tubeworm-associated samples, diversity is calculated individually on all extracts and then with size fractions and method replicates combined. Wilcoxon Rank Sum tests are applied to test for significant differences in diversity levels.
3.1 Macrofauna
MacroDiv <- diversity(macro9, index = "invsimpson", MARGIN = 1)
MacroDiv## EMw3 ECw11 EMw1 EMw4 EMw8 MVw12 EMw5 ECw10
## 2.357427 2.626470 3.031563 3.046310 3.043524 2.297812 1.266767 1.940417
## MVw13
## 9.226349Wilcox test between highT and lowT diversity levels.
MAChighTdiv <- c(3.03, 3.04, 3.04, 2.3, 1.94)
MAClowTdiv <- c(2.36, 2.62, 1.27, 9.23)
wilcox.test(MAClowTdiv, MAChighTdiv)##
## Wilcoxon rank sum test with continuity correction
##
## data: MAClowTdiv and MAChighTdiv
## W = 9, p-value = 0.9021
## alternative hypothesis: true location shift is not equal to 0Macrofaunal diversity IS NOT significantly different between highT and lowT grabs.
3.2 Meiofauna
MeioDiv <- diversity(meio9, index = "invsimpson", MARGIN = 1)
MeioDiv## EMw3 ECw11 EMw1 EMw4 EMw8 MVw12 EMw5 ECw10
## 4.784573 4.867552 1.242057 1.078020 1.012415 1.000000 3.020960 1.380738
## MVw13
## 1.543465Wilcox test between highT and lowT diversity levels.
MEhighTdiv <- c(1.24, 1.07, 1.01, 1, 1.38)
MElowTdiv <- c(4.78, 4.86, 3.02, 1.54)
wilcox.test(MEhighTdiv, MElowTdiv)##
## Wilcoxon rank sum exact test
##
## data: MEhighTdiv and MElowTdiv
## W = 0, p-value = 0.01587
## alternative hypothesis: true location shift is not equal to 0Meiofaunal diversity IS significantly different between highT and lowT grabs.
3.3 Microbes (40 extracts)
For individual DNA extracts, sample naming conventions for tubeworm associated extracts are as follows:
Sample name followed by microbial community removal method (bk= agitation in bucket of artificial seawater, d= direct removal of attached biofilm) and size fraction (m= micro 20-64µm, pn= pico/nano 0.2-20µm).
EXAMPLE: EMw3bkm –> sample EMw3, bucket agitation method, micro size fraction.
3.3.1 Bacteria
Bact40nonsingDiv <- diversity(Bact40, index = "invsimpson", MARGIN = 1)
Bact40nonsingDiv## EMw3bkm EMw3bkpn EMf3 EMw6bkm EMw6bkpn EMw6dm EMw6dpn
## 13.877373 132.019034 3.356216 190.791269 129.672086 96.763474 180.551378
## ECw11dm ECw11dpn ECf11 EMw1bkm EMw1bkpn EMw1dm EMw1dpn
## 18.147660 12.664641 13.794946 69.317711 67.080719 13.660772 22.860723
## EMf1 EMw4dm EMw4dpn EMw7dm EMw7dpn EMf7 EMw8dm
## 6.725429 5.003208 5.384283 30.205731 16.667821 8.739932 26.253335
## EMw8dpn MVw12bkm MVw12bkpn MVf12 EMw2dm EMw2dpn EMf2
## 30.179900 11.177898 13.907539 6.783713 114.013041 61.603672 6.001158
## EMw5bkm EMw5bkpn ECw9bkm ECw9bkpn ECf9 ECw10bkm ECw10bkpn
## 108.581991 13.380956 38.195485 125.157585 10.439104 10.566329 15.878173
## MVw13bkm MVw13bkpn bkCB bkEM bkMV
## 30.641508 39.521918 9.444129 13.949291 16.312045BhighTdiv <- c(11.18, 10.57, 69.32, 30.21, 13.66, 5, 26.25, 13.91, 15.88,
67.08, 16.67, 22.86, 5.38, 30.18)
BlowTdiv <- c(190.79, 13.88, 108.58, 38.2, 30.64, 96.76, 18.15, 114.01,
129.67, 132.02, 13.38, 125.16, 39.52, 180.55, 12.66, 61.6)
BfluidsDiv <- c(9.44, 13.95, 16.31, 3.36, 13.79, 6.73, 8.74, 6.78, 6, 10.44)
wilcox.test(BhighTdiv, BlowTdiv)##
## Wilcoxon rank sum exact test
##
## data: BhighTdiv and BlowTdiv
## W = 43, p-value = 0.003321
## alternative hypothesis: true location shift is not equal to 0wilcox.test(BfluidsDiv, BhighTdiv)##
## Wilcoxon rank sum exact test
##
## data: BfluidsDiv and BhighTdiv
## W = 30, p-value = 0.01855
## alternative hypothesis: true location shift is not equal to 0wilcox.test(BfluidsDiv, BlowTdiv)##
## Wilcoxon rank sum exact test
##
## data: BfluidsDiv and BlowTdiv
## W = 8, p-value = 2.523e-05
## alternative hypothesis: true location shift is not equal to 0Bacterial diversity IS signficantly different between highT and lowT grabs and between fluids and both highT and lowT grabs.
3.3.2 Eukarya
Euk40nonsingDiv <- diversity(Euk40, index = "invsimpson", MARGIN = 1)
Euk40nonsingDiv## EMw3bkm EMw3bkpn EMf3 EMw6bkm EMw6bkpn EMw6dm EMw6dpn ECw11dm
## 24.935378 7.484512 20.319558 11.688583 2.682649 2.561250 3.220499 9.205386
## ECw11dpn ECf11 EMw1bkm EMw1bkpn EMw1dm EMw1dpn EMf1 EMw4dm
## 3.037183 7.435008 8.593577 2.427560 1.325924 1.225294 5.474317 4.348228
## EMw4dpn EMw7dm EMw7dpn EMf7 EMw8dm EMw8dpn MVw12bkpn MVw12bm
## 3.277309 2.936878 3.693987 14.214942 2.701855 5.718646 11.698233 5.852115
## MVf12 EMw2dm EMw2dpn EMf2 EMw5bkm EMw5bkpn ECw9bkm ECw9bkpn
## 8.343386 8.200491 8.587284 16.754450 27.093450 40.638687 8.334630 17.235464
## ECf9 ECw10bkm ECw10bkpn MVw13bkm MVw13bkpn bkCB bkEM bkMV
## 1.751818 4.597617 13.682620 6.739444 8.104348 2.013199 1.829072 33.575884EhighTdiv <- c(11.7, 4.6, 2.94, 8.59, 1.33, 4.35, 2.7, 5.85, 13.68, 3.69,
2.43, 1.23, 3.28, 5.72)
ElowTdiv <- c(11.69, 2.56, 24.94, 9.21, 8.2, 27.09, 8.33, 6.74, 2.68, 3.22,
7.48, 3.04, 8.59, 40.64, 17.24, 8.1)
EfluidsDiv <- c(2.01, 1.83, 33.58, 20.32, 7.44, 5.47, 14.21, 8.34, 16.75,
1.75)
wilcox.test(EhighTdiv, ElowTdiv)##
## Wilcoxon rank sum test with continuity correction
##
## data: EhighTdiv and ElowTdiv
## W = 61.5, p-value = 0.03764
## alternative hypothesis: true location shift is not equal to 0wilcox.test(EfluidsDiv, EhighTdiv)##
## Wilcoxon rank sum exact test
##
## data: EfluidsDiv and EhighTdiv
## W = 93, p-value = 0.1915
## alternative hypothesis: true location shift is not equal to 0wilcox.test(ElowTdiv, EfluidsDiv)##
## Wilcoxon rank sum exact test
##
## data: ElowTdiv and EfluidsDiv
## W = 90, p-value = 0.6227
## alternative hypothesis: true location shift is not equal to 0Eukaryal diversity IS signficantly different between highT and lowT grabs, but not between fluids and either highT or lowT grabs.
3.3.3 Archaea
Arch40nonsingDiv <- diversity(Arch38, index = "invsimpson", MARGIN = 1)
Arch40nonsingDiv## EMw3bkm EMw3bkpn EMf3pn EMw6bkm EMw6bkpn EMw6dm EMw6dpn ECf11pn
## 1.020409 1.034489 4.710203 1.042860 1.032802 1.182467 1.012818 3.541635
## EMw1bkm EMw1bkpn EMw1dm EMw1dpn EMf1pn EMw4dm EMw4dpn EMw7dm
## 1.121875 1.619707 6.671790 6.443250 3.614836 2.096764 6.405379 6.050011
## EMw7dpn EMf7pn EMw8dm EMw8dpn MVw12bkm MVw12bkpn MVf12pn EMw2dm
## 7.688165 3.067253 5.994080 4.717821 3.940932 7.645768 4.493723 2.341680
## EMw2dpn EMf2pn EMw5bkm EMw5bkpn ECw9bkm ECw9bkpn ECf9pn ECw10bkm
## 1.030438 4.646234 1.018484 1.014421 4.146493 2.351466 4.397208 8.000000
## ECw10bkpn MVw13bkm MVw13bkpn bkECpn bkEMpn bkMVpn
## 9.403731 2.860098 2.608518 3.242648 2.936804 5.076730AhighTdiv <- c(3.94, 8, 6.05, 1.12, 6.67, 2.1, 5.99, 7.65, 9.4, 7.69, 1.62,
6.44, 6.41, 4.72)
AlowTdiv <- c(1.04, 1.18, 1.02, 2.34, 1.02, 4.15, 2.86, 1.03, 1.01, 1.03,
1.03, 1.01, 2.35, 2.61)
AfluidsDiv <- c(3.24, 2.94, 5.08, 4.71, 3.54, 3.61, 3.07, 4.49, 4.65, 4.4)
wilcox.test(AlowTdiv, AhighTdiv)##
## Wilcoxon rank sum test with continuity correction
##
## data: AlowTdiv and AhighTdiv
## W = 17, p-value = 0.0002141
## alternative hypothesis: true location shift is not equal to 0wilcox.test(AfluidsDiv, AhighTdiv)##
## Wilcoxon rank sum exact test
##
## data: AfluidsDiv and AhighTdiv
## W = 36, p-value = 0.0484
## alternative hypothesis: true location shift is not equal to 0wilcox.test(AlowTdiv, AfluidsDiv)##
## Wilcoxon rank sum test with continuity correction
##
## data: AlowTdiv and AfluidsDiv
## W = 5, p-value = 0.0001558
## alternative hypothesis: true location shift is not equal to 0Archaeal diversity IS signficantly different between highT and lowT grabs and between fluids and both highT and lowT grabs.
3.4 Diversity boxplots
(Plotting for all 40 extracts)
# Combine diversity values into one dataframe
InvSimpson <- c(BhighTdiv, BlowTdiv, BfluidsDiv, EhighTdiv, ElowTdiv, EfluidsDiv,
AhighTdiv, AlowTdiv, AfluidsDiv, MAChighTdiv, MAClowTdiv, MEhighTdiv,
MElowTdiv)
domain <- c(rep("Bacteria", 40), rep("Eukarya", 40), rep("Archaea", 38),
rep("Macrofauna", 9), rep("Meiofauna", 9))
Sample_type <- c(rep(c(rep("highT", 14), rep("lowT", 16), rep("fluids",
10)), 2), rep("highT", 14), rep("lowT", 14), rep("fluids", 10), rep(c(rep("highT",
5), rep("lowT", 4)), 2))
Div4Plots <- data.frame(InvSimpson, Sample_type, domain)
# Plot them
ggplot(Div4Plots, aes(domain, InvSimpson, fill = Sample_type)) + geom_boxplot(varwidth = FALSE) +
scale_fill_discrete(breaks = c("fluids", "highT", "lowT"), labels = c("Fluids",
"HighT", "LowT"), name = "Sample Type") + labs(x = NULL, y = "Diversity (Inverse Simpson)") +
facet_wrap(. ~ domain, scales = "free") + theme(axis.text.x = element_blank(),
axis.title.x = element_blank(), axis.text.y = element_text(size = 12),
axis.title.y = element_text(size = 14))
3.5 Microbes (23 samples)
For these calculations, size fractions and method replicates are combined.
3.5.1 Bacteria
NOTE: These diversity calculations are based on subsampled datasets (using the smallest read count by domain).
rowSums(codaBact23nonsingleton)## EMf3 ECf11 EMf1 EMf7 MVf12 EMf2 ECf9 bkEC bkEM bkMV EMw3
## 34535 20711 28350 30036 18755 34250 24775 21278 34180 26397 44531
## EMw6 ECw11 EMw1 EMw4 EMw7 EMw8 MVw12 EMw2 EMw5 ECw9 ECw10
## 83006 32677 112883 58955 60625 64972 45872 66446 30031 49338 51516
## MVw13
## 24247subBact23 <- as.data.frame(rrarefy(codaBact23nonsingleton, sample = 18755))
subBact23Div <- diversity(subBact23, index = "invsimpson", MARGIN = 1)
subBact23Div## EMf3 ECf11 EMf1 EMf7 MVf12 EMf2 ECf9
## 3.380745 13.753857 6.662416 8.823526 6.783713 5.940150 10.509222
## bkEC bkEM bkMV EMw3 EMw6 ECw11 EMw1
## 9.462509 14.124186 16.546602 36.344101 217.909681 16.925666 47.139049
## EMw4 EMw7 EMw8 MVw12 EMw2 EMw5 ECw9
## 5.394489 25.912545 27.535837 13.170838 94.589946 95.946261 73.961936
## ECw10 MVw13
## 12.918118 37.5450773.5.2 Eukarya
subEuk23 <- as.data.frame(rrarefy(codaEuk23nonsingleton, sample = 13172))
subEuk23Div <- diversity(subEuk23, index = "invsimpson", MARGIN = 1)
subEuk23Div## EMf3 ECf11 EMf1 EMf7 MVf12 EMf2 ECf9 bkEC
## 21.132670 7.371347 5.328299 14.430674 8.329511 17.015357 1.772109 1.984034
## bkEM bkMV EMw3 EMw6 ECw11 EMw1 EMw4 EMw7
## 1.863651 34.021920 8.507661 5.353393 8.783340 4.391468 6.282467 4.685926
## EMw8 MVw12 EMw2 EMw5 ECw9 ECw10 MVw13
## 6.034257 8.896173 10.016093 42.089294 15.393216 10.769617 7.6045503.5.3 Archaea
subArch22 <- as.data.frame(rrarefy(codaArch22nonsingleton, sample = 25064))
subArch22Div <- diversity(subArch22, index = "invsimpson", MARGIN = 1)
subArch22Div## EMf3 ECf11 EMf1 EMf7 MVf12 EMf2 ECf9 bkEC
## 4.711096 3.512507 3.605891 3.024852 4.461340 4.650887 4.486183 3.278287
## bkEM bkMV EMw3 EMw6 EMw1 EMw4 EMw7 EMw8
## 2.907239 5.035000 1.028663 1.025469 1.805728 5.709097 6.991795 4.897325
## MVw12 EMw2 EMw5 ECw9 ECw10 MVw13
## 5.407972 1.031313 1.014593 2.875343 9.411689 2.7316274 Comparing sample compositions
4.1 Compositional data approach
Data from high-throughput sequencing platforms, which have an arbitrary sum imposed by the instrument, contain only information on ratios between taxa and are therefore considered compositional in nature. Data analysis followed the compositional approach recommended by Gloor and Reid (2016) and Quinn et al (2019). Zero counts were replaced by an imputed value using the count zero multiplicative method in the zCompositions R package (Martín-Fernández et al 2011, Palarea-Albaladejo and Martin-Fernandez 2015), and data were then transformed by centered log-ratio (Aitchison 1986). Faunal data followed the same approach because full assemblages were not collected.
4.2 Macrofauna
Nonmetric Multidimensional Scaling to visualize the relationships between samples based on species abundances.
# Package 'zCompositions' for imputation of zeros. Output='p-counts'
# returns pseudo-counts (preserves ratios in the original data after
# conversion of zeros to non-zeros).
mac9_cmultRepl <- cmultRepl(macro9, method = "CZM", output = "p-counts")
# Centered log-ratio transformation
clr_mac9 <- apply(mac9_cmultRepl, 2, function(x) {
log(x) - mean(log(x))
})
# Package 'compositions' to calculate Aitchison distance matrix metaMDS
# for Nonmetric Multidimensional Scaling
clr_mac9.dist <- dist(clr_mac9)
clr_mac9.NMDS <- metaMDS(clr_mac9.dist, k = 2, try = 100)
# Trim the metadata file down to only the samples with all three size
# classes characterized and re-order the sample rows to match the
# macro9 file.
metadata9 <- metadata23[c(1, 3, 4, 5, 8, 10, 11, 12, 13), ]
metadata9 <- metadata9[match(rownames(macro9), rownames(metadata9)), ]
# Make a dataframe with NMDS points and associated metadata
clr_mac9.NMDS.df <- data.frame(x = clr_mac9.NMDS$points[, 1], y = clr_mac9.NMDS$points[,
2], T_category = as.factor(metadata9[, 5]), Vent_field = as.factor(metadata9[,
3]), Temperature = metadata9[, 4], Substratum = as.factor(metadata9[,
6]))
# Plot it
ggplot(data = clr_mac9.NMDS.df, aes(x = x, y = y, color = Temperature,
shape = Substratum)) + geom_point(size = 4) + geom_text_repel(aes(label = rownames(metadata9)),
hjust = 0.5, vjust = 1.5, size = 4) + ggtitle("Macrofauna NMDS") +
scale_color_gradient(low = "blue", high = "red")
Envfit to see which variables (Vent_field + Substratum + Temperature + T_category) explain the NMDS plot.
Temperature is a signficant predictor, followed by substratum.
clr_mac9_NMDSfit <- envfit(clr_mac9.NMDS ~ ., metadata9)
clr_mac9_NMDSfit##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## Temperature 0.99151 -0.13002 0.927 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## NMDS1 NMDS2
## SourceTubeworms 0.0000 0.0000
## Sample_typebk -0.7698 -0.0033
## Sample_typed 1.5396 0.0066
## Vent_fieldClam_Bed -0.7717 -0.0030
## Vent_fieldMain_Endeavour 0.6161 0.0021
## Vent_fieldMiddle_Valley -0.7686 -0.0021
## T_categoryhigh 4.7975 -2.3862
## T_categorylow -5.9969 2.9827
## Substratumbasalt -7.7221 -0.0325
## Substratumsulfide 2.2063 0.0093
##
## Goodness of fit:
## r2 Pr(>r)
## Source 0.0000 1.000
## Sample_type 0.0171 1.000
## Vent_field 0.0069 0.988
## T_category 0.5192 0.012 *
## Substratum 0.2465 0.026 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999Hierarchical cluster analysis.
clr_mac9.CLUST <- hclust(clr_mac9.dist, method = "average", member = NULL)
plot(clr_mac9.CLUST, hang = -1, main = "Macrofauna")
ANOSIM tests on the clr-transformed data to assess differences based on highT vs lowT, substratum and vent field.
mac9ANOSIMtemp <- with(metadata9, anosim(clr_mac9.dist, T_category))
summary(mac9ANOSIMtemp)##
## Call:
## anosim(x = clr_mac9.dist, grouping = T_category)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: 0.7812
## Significance: 0.012
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.269 0.325 0.381 0.781
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 13 17.75 25.5 31.25 36 20
## high 1 3.25 6.0 10.75 23 10
## low 6 8.25 14.5 21.50 28 6mac9ANOSIMsubstr <- with(metadata9, anosim(clr_mac9.dist, Substratum))
summary(mac9ANOSIMsubstr)##
## Call:
## anosim(x = clr_mac9.dist, grouping = Substratum)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: 0.3961
## Significance: 0.107
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.396 0.474 0.656 0.656
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 6 17.25 24.5 29.75 35 14
## basalt 8 8.00 8.0 8.00 8 1
## sulfide 1 7.00 14.0 23.00 36 21mac9ANOSIMfield <- with(metadata9, anosim(clr_mac9.dist, Vent_field))
summary(mac9ANOSIMfield)##
## Call:
## anosim(x = clr_mac9.dist, grouping = Vent_field)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: 0.2847
## Significance: 0.091
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.278 0.368 0.458 0.646
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 4 10.75 21.5 28.25 35 24
## Clam_Bed 32 32.00 32.0 32.00 32 1
## Main_Endeavour 1 3.75 13.5 15.75 25 10
## Middle_Valley 36 36.00 36.0 36.00 36 1RESULT: Temperature category is a significant factor in macrofaunal diversity variation, substratum and vent field are not.
4.3 Meiofauna
Envfit. Temperature is the only signficant predictor.
##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## Temperature 0.99407 0.10876 0.9072 0.003 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## NMDS1 NMDS2
## SourceTubeworms 0.0000 0.0000
## Sample_typebk -1.0118 1.4433
## Sample_typed 2.0236 -2.8866
## Vent_fieldClam_Bed -3.1965 -3.2965
## Vent_fieldMain_Endeavour -0.1757 -3.0498
## Vent_fieldMiddle_Valley 3.6357 10.9209
## T_categoryhigh 12.5709 -2.0716
## T_categorylow -15.7136 2.5895
## Substratumbasalt -19.1619 -4.5150
## Substratumsulfide 5.4748 1.2900
##
## Goodness of fit:
## r2 Pr(>r)
## Source 0.0000 1.000
## Sample_type 0.0218 0.700
## Vent_field 0.1376 0.738
## T_category 0.7105 0.008 **
## Substratum 0.3877 0.054 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
ANOSIM tests on the clr-transformed data (highT vs lowT, substratum and vent field).
##
## Call:
## anosim(x = clr_meio9.dist, grouping = T_category)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: 0.9562
## Significance: 0.004
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.344 0.425 0.444 0.669
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 15 21.50 26.5 31.25 36 20
## high 1 3.25 5.5 10.00 13 10
## low 8 9.25 12.0 16.25 21 6##
## Call:
## anosim(x = clr_meio9.dist, grouping = Substratum)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: 0.5325
## Significance: 0.06
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.175 0.532 0.695 0.695
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 8 19.5 25.5 31.5 36 14
## basalt 10 10.0 10.0 10.0 10 1
## sulfide 1 6.0 15.0 22.0 33 21##
## Call:
## anosim(x = clr_meio9.dist, grouping = Vent_field)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: -0.1806
## Significance: 0.862
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.285 0.396 0.467 0.507
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 1 9.75 17.5 23.25 36 24
## Clam_Bed 30 30.00 30.0 30.00 30 1
## Main_Endeavour 3 8.75 25.5 28.75 34 10
## Middle_Valley 15 15.00 15.0 15.00 15 1RESULT: Temperature category is a significant factor in meiofaunal diversity variation, substratum and vent field are not.
4.4 Microbes (9 fully characterized samples)
4.4.1 Bacteria
codaBact9 <- codaBact23nonsingleton[c(11, 13, 14, 15, 17, 18, 20, 22, 23),
]
# remove any OTU with only 1 read. Drops from 40,389 to 18,570 OTUs.
codaBact9nonsingleton <- codaBact9[, which(colSums(codaBact9) > 1)]
Envfit. Temperature is the only signficant predictor.
##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## Temperature 0.67513 0.73770 0.8323 0.005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## NMDS1 NMDS2
## SourceTubeworms 0.0000 0.0000
## Sample_typebk -9.1711 -10.3184
## Sample_typed 18.3422 20.6368
## Vent_fieldClam_Bed -3.0742 -3.6403
## Vent_fieldMain_Endeavour 16.4566 -17.8509
## Vent_fieldMiddle_Valley -38.0674 48.2674
## T_categoryhigh 38.6052 21.1261
## T_categorylow -48.2565 -26.4076
## Substratumbasalt -46.6324 -58.8916
## Substratumsulfide 13.3235 16.8262
##
## Goodness of fit:
## r2 Pr(>r)
## Source 0.0000 1.000
## Sample_type 0.0569 0.734
## Vent_field 0.1749 0.680
## T_category 0.3612 0.017 *
## Substratum 0.2405 0.169
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
ANOSIM tests on the clr-transformed data (Temp category and substratum).
##
## Call:
## anosim(x = clr_bact9.dist, grouping = Substratum)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: 0.1494
## Significance: 0.342
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.448 0.513 0.838 0.838
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 4 11.75 21.5 29.75 34 14
## basalt 22 22.00 22.0 22.00 22 1
## sulfide 1 9.00 17.0 25.00 36 21##
## Call:
## anosim(x = clr_bact9.dist, grouping = T_category)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: 0.325
## Significance: 0.03
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.175 0.275 0.325 0.381
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 5 14.00 19.5 30.25 35 20
## high 1 5.25 17.0 24.75 36 10
## low 4 7.50 11.5 20.00 23 6RESULT: Temperature category is a weak factor in bacterial diversity variation, substratum is not.
4.4.2 Eukarya
Envfit. No signficant predictors.
##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## Temperature -0.05364 0.99856 0.5827 0.093 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## NMDS1 NMDS2
## SourceTubeworms 0.0000 0.0000
## Sample_typebk 2.2839 -7.0991
## Sample_typed -4.5677 14.1982
## Vent_fieldClam_Bed 9.6160 9.8077
## Vent_fieldMain_Endeavour 9.3395 -13.1135
## Vent_fieldMiddle_Valley -32.9648 22.9760
## T_categoryhigh 1.4697 21.6364
## T_categorylow -1.8372 -27.0455
## Substratumbasalt 19.4432 -33.6297
## Substratumsulfide -5.5552 9.6085
##
## Goodness of fit:
## r2 Pr(>r)
## Source 0.0000 1.000
## Sample_type 0.0257 0.859
## Vent_field 0.1259 0.818
## T_category 0.1359 0.547
## Substratum 0.0997 0.597
## Permutation: free
## Number of permutations: 999
ANOSIM tests on the clr-transformed data (Temp category and substratum).
##
## Call:
## anosim(x = clr_euk9.dist, grouping = Substratum)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: 0.2922
## Significance: 0.265
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.558 0.695 0.825 0.825
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 11 17.25 21.5 26.5 33 14
## basalt 26 26.00 26.0 26.0 26 1
## sulfide 1 6.00 12.0 30.0 36 21##
## Call:
## anosim(x = clr_euk9.dist, grouping = T_category)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: 0.05625
## Significance: 0.291
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.175 0.225 0.256 0.312
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 4 12.75 18.5 24.25 36 20
## high 1 3.75 8.0 30.75 34 10
## low 11 17.00 24.5 27.50 29 6RESULT: Neither temperature nor substratum affect variation in microeukaryote diveristy.
4.4.3 Archaea
Sample ECw11 produced insuffient archaeal reads so there are only 8 samples included here.
Envfit. No signficant predictors.
##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## Temperature -7.3631e-05 -1.0000e+00 0.3679 0.303
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## NMDS1 NMDS2
## SourceTubeworms 0.0000 0.0000
## Sample_typebk 4.1121 0.0005
## Sample_typed -12.3363 -0.0015
## Vent_fieldClam_Bed -12.3328 -0.0156
## Vent_fieldMain_Endeavour -12.3347 0.0012
## Vent_fieldMiddle_Valley 37.0032 0.0047
## T_categoryhigh -12.3359 -0.0038
## T_categorylow 20.5598 0.0063
## Substratumbasalt -12.3271 0.0094
## Substratumsulfide 1.7610 -0.0013
##
## Goodness of fit:
## r2 Pr(>r)
## Source 0.0000 1.000
## Sample_type 0.0476 0.588
## Vent_field 0.4286 0.266
## T_category 0.2382 0.063 .
## Substratum 0.0204 0.890
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
ANOSIM tests on the clr-transformed data (Temp category and substratum).
##
## Call:
## anosim(x = clr_arch8.dist, grouping = Substratum)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: -0.483
## Significance: 1
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 1 1 1 1
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 1 3 8 14.5 22 7
## basalt NA NA NA NA NA 0
## sulfide 3 10 16 23.0 28 21##
## Call:
## anosim(x = clr_arch8.dist, grouping = T_category)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: -0.01538
## Significance: 0.476
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.282 0.364 0.405 0.569
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 2 5.50 12.0 24.50 28 15
## high 7 10.75 14.5 18.25 21 10
## low 1 11.50 22.0 22.50 23 3RESULT: Neither temperature or substratum have an effect on variation in archaeal diversity.
4.5 Microbes (all 23 samples)
4.5.1 Three domains, QPCR-balanced
Rel
Envfit. Temperature is the only significant predictor.
##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## Temperature 0.15956 -0.98719 0.0717 0.489
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## NMDS1 NMDS2
## SourceBackground -7.4027 29.4711
## SourceDiffuse 5.6560 3.7974
## SourceTubeworms -1.3373 -8.8458
## Vent_fieldClam_Bed 26.0750 -8.7773
## Vent_fieldMain_Endeavour -34.6422 -1.8391
## Vent_fieldMiddle_Valley 73.4745 19.1430
## Substratumbasalt 20.8456 -3.7023
## Substratumnone -7.4027 29.4711
## Substratumsulfide -5.4680 -4.6601
##
## Goodness of fit:
## r2 Pr(>r)
## Source 0.0282 0.902
## Vent_field 0.2953 0.014 *
## Substratum 0.0395 0.819
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
ANOSIM tests (Vent field, substratum, sample type, temperature category)
##
## Call:
## anosim(x = clr_balMic.dist, grouping = Vent_field)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: 0.4485
## Significance: 0.002
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.169 0.225 0.263 0.319
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 1 99.75 155.0 209.75 253 154
## Clam_Bed 3 83.00 175.0 225.50 246 15
## Main_Endeavour 2 26.50 68.5 111.50 184 78
## Middle_Valley 120 177.50 180.0 187.75 204 6##
## Call:
## anosim(x = clr_balMic.dist, grouping = Source)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: 0.2771
## Significance: 0.021
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.144 0.192 0.251 0.314
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 13 93.50 136 186.0 252 151
## Background 36 43.00 50 58.0 66 3
## Diffuse 4 83.00 174 222.0 253 21
## Tubeworms 1 25.25 53 187.5 236 78##
## Call:
## anosim(x = clr_balMic.dist, grouping = Substratum)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: 0.1286
## Significance: 0.155
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.178 0.242 0.302 0.389
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 2 74.50 138.0 199.00 252 135
## basalt 15 189.75 230.5 246.25 253 10
## none 36 43.00 50.0 58.00 66 3
## sulfide 1 52.00 114.0 169.00 241 105##
## Call:
## anosim(x = clr_balMic.dist, grouping = T_category)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: -0.05103
## Significance: 0.584
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.185 0.276 0.347 0.407
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 1 59.25 120.5 185.00 252 102
## high 6 18.50 26.0 207.50 236 15
## low 4 75.75 135.0 190.25 253 136RESULT: Variation in microbial community diversity is not tied to any of the measured variables.
4.5.2 Bacteria
# Imputation of zeros.
codaBact23_cmultRepl <- cmultRepl(codaBact23nonsingleton, method = "CZM",
output = "p-counts")
# Need metadata for the next part. Must be in the same order as the
# samples.
metadata23 <- metadata23[match(rownames(codaBact23nonsingleton), rownames(metadata23)),
]
# Centered log-ratio transformation
clr_codaBact23 <- apply(codaBact23_cmultRepl, 2, function(x) {
log(x) - mean(log(x))
})
# Calculate Aitchison distance matrix
clr_codaBact23.dist <- dist(clr_codaBact23)
# NMDS
clr_codaBact23.NMDS <- metaMDS(clr_codaBact23.dist, k = 2)
clr_codaBact23.NMDS.df <- data.frame(x = clr_codaBact23.NMDS$points[, 1],
y = clr_codaBact23.NMDS$points[, 2], Source = as.factor(metadata23[,
1]), Vent_field = as.factor(metadata23[, 3]), Temperature = metadata23[,
4], T_category = as.factor(metadata23[, 5]), Substratum = as.factor(metadata23[,
6]))
ggplot(data = clr_codaBact23.NMDS.df, aes(x = x, y = y, color = Temperature,
shape = Source)) + geom_point(size = 4) + geom_text(aes(label = rownames(metadata23)),
hjust = 0.5, vjust = 1.5, size = 4) + ggtitle("Bacteria (23)") + scale_color_gradient(low = "blue",
high = "red")
Envfit. Temperature is the one and only significant predictor.
##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## Temperature -0.83099 0.55628 0.485 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## NMDS1 NMDS2
## Vent_fieldClam_Bed -4.6834 -1.0941
## Vent_fieldMain_Endeavour 5.9497 7.3926
## Vent_fieldMiddle_Valley -12.3115 -22.3849
## Substratumbasalt 37.8105 13.7544
## Substratumnone -7.2824 -24.7856
## Substratumsulfide -11.1470 0.3723
## SourceBackground -7.2824 -24.7856
## SourceDiffuse -6.7661 -20.9488
## SourceTubeworms 5.3238 16.9999
##
## Goodness of fit:
## r2 Pr(>r)
## Vent_field 0.0424 0.806
## Substratum 0.1294 0.221
## Source 0.1030 0.358
## Permutation: free
## Number of permutations: 999ANOSIM tests (Vent field, substratum, sample type, temperature category)
##
## Call:
## anosim(x = clr_codaBact23.dist, grouping = Vent_field)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: -0.1619
## Significance: 0.937
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.171 0.222 0.270 0.332
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 3 58.25 117.5 177.75 247 154
## Clam_Bed 2 65.50 85.0 117.00 167 15
## Main_Endeavour 1 102.25 172.0 216.75 253 78
## Middle_Valley 31 43.00 60.5 85.50 109 6##
## Call:
## anosim(x = clr_codaBact23.dist, grouping = Substratum)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: -0.04545
## Significance: 0.563
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.206 0.268 0.324 0.423
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 1 55.50 120 189.00 253 135
## basalt 32 86.25 166 195.75 241 10
## none 3 5.50 8 11.00 14 3
## sulfide 4 83.00 131 189.00 252 105##
## Call:
## anosim(x = clr_codaBact23.dist, grouping = Source)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: -0.07129
## Significance: 0.682
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.154 0.206 0.267 0.323
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 2 68.5 115.0 170.5 247 151
## Background 3 5.5 8.0 11.0 14 3
## Diffuse 1 9.0 18.0 32.0 53 21
## Tubeworms 40 123.5 184.5 209.5 253 78##
## Call:
## anosim(x = clr_codaBact23.dist, grouping = T_category)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: 0.4051
## Significance: 0.011
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.210 0.281 0.333 0.398
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 70 110.75 148.5 207.25 253 102
## high 40 54.00 119.0 186.50 233 15
## low 1 34.75 79.0 174.25 247 136RESULT: Variation in bacterial diversity is tied to temperature, but not to vent field, substratum or sample type.
4.5.3 Eukarya
Envfit. Sample type is a weak predictor.
##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## Temperature 0.12405 0.99228 0.0231 0.803
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## NMDS1 NMDS2
## Vent_fieldClam_Bed 6.8767 4.6585
## Vent_fieldMain_Endeavour 8.3035 -5.6479
## Vent_fieldMiddle_Valley -37.3013 11.3678
## Substratumbasalt 20.3966 -4.5481
## Substratumnone -20.2410 -6.4232
## Substratumsulfide -2.7507 2.8007
## SourceBackground -20.2410 -6.4232
## SourceDiffuse -25.4307 -32.3771
## SourceTubeworms 18.3645 18.9161
##
## Goodness of fit:
## r2 Pr(>r)
## Vent_field 0.0872 0.446
## Substratum 0.0421 0.771
## Source 0.2486 0.022 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999ANOSIM tests (Vent field, substratum, sample type, temperature category)
##
## Call:
## anosim(x = clr_codaEuk23.dist, grouping = Vent_field)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: -0.04486
## Significance: 0.617
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.156 0.223 0.254 0.309
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 1 62.25 123.0 192.75 253 154
## Clam_Bed 5 39.00 86.0 168.50 222 15
## Main_Endeavour 2 75.50 145.0 194.75 252 78
## Middle_Valley 95 100.50 108.5 136.00 199 6##
## Call:
## anosim(x = clr_codaEuk23.dist, grouping = Substratum)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: 0.03239
## Significance: 0.379
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.197 0.247 0.314 0.374
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 4 71.0 124 192.0 253 135
## basalt 14 154.5 192 227.5 252 10
## none 1 26.5 52 57.5 63 3
## sulfide 2 59.0 127 185.0 247 105##
## Call:
## anosim(x = clr_codaEuk23.dist, grouping = Source)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: 0.1153
## Significance: 0.15
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.149 0.208 0.245 0.324
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 4 77.50 131 192.50 252 151
## Background 1 26.50 52 57.50 63 3
## Diffuse 8 18.00 38 65.00 166 21
## Tubeworms 2 79.25 148 200.75 253 78##
## Call:
## anosim(x = clr_codaEuk23.dist, grouping = T_category)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: -0.1266
## Significance: 0.818
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.185 0.249 0.289 0.344
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 12 67.25 102.5 159.50 247 102
## high 2 8.50 27.0 167.50 221 15
## low 1 81.75 147.5 198.25 253 136RESULT: Variation in microeukaryote diversity is not tied to any of the measured variables.
4.5.4 Archaea
Envfit. Temperature and sample type are significant predictors.
##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## Temperature 0.68773 0.72597 0.603 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## NMDS1 NMDS2
## Vent_fieldClam_Bed 3.5115 -8.8294
## Vent_fieldMain_Endeavour -4.3847 7.6991
## Vent_fieldMiddle_Valley 9.8609 -13.9853
## Substratumbasalt -11.0887 -8.9927
## Substratumnone -18.6048 0.7761
## Substratumsulfide 6.6780 2.2428
## SourceBackground -18.6048 0.7761
## SourceDiffuse -30.7007 -7.6776
## SourceTubeworms 22.5600 4.2846
##
## Goodness of fit:
## r2 Pr(>r)
## Vent_field 0.0820 0.487
## Substratum 0.0806 0.512
## Source 0.4460 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999ANOSIM tests (Vent field, substratum, sample type, temperature category)
##
## Call:
## anosim(x = clr_codaArch22.dist, grouping = Vent_field)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: 0.06507
## Significance: 0.297
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.170 0.213 0.258 0.312
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 7 58.00 126.0 179.00 231 137
## Clam_Bed 22 53.75 156.5 197.25 228 10
## Main_Endeavour 1 60.75 103.5 156.75 226 78
## Middle_Valley 10 69.75 121.0 172.25 178 6##
## Call:
## anosim(x = clr_codaArch22.dist, grouping = Source)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: 0.2482
## Significance: 0.02
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.157 0.200 0.232 0.279
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 9 62.00 136.0 186.00 231 141
## Background 7 7.50 8.0 10.50 13 3
## Diffuse 14 52.00 79.0 99.00 167 21
## Tubeworms 1 63.75 113.5 152.75 218 66##
## Call:
## anosim(x = clr_codaArch22.dist, grouping = Substratum)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: -0.2803
## Significance: 0.985
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.198 0.269 0.326 0.360
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 1 44.0 90 155.0 227 117
## basalt 4 66.0 81 93.0 111 6
## none 7 7.5 8 10.5 13 3
## sulfide 5 94.0 146 194.0 231 105##
## Call:
## anosim(x = clr_codaArch22.dist, grouping = T_category)
## Dissimilarity: euclidean
##
## ANOSIM statistic R: 0.4778
## Significance: 0.001
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.167 0.228 0.262 0.340
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 23 112.75 157.0 192.25 231 96
## high 45 103.50 134.0 156.50 218 15
## low 1 34.75 74.5 138.50 230 120RESULT: Variation in archaeal diversity is tied to temperature and sample type.
5 Relative Abundance plots
Convert counts to relatvie abundances for stacked bar plots.
5.1 Macrofauna
Group “Others” include taxa that are always less than 2% relative abundance.
# Convert count data to relative abundances (it must be as a matrix)
macro9matrix <- as.matrix(macro9)
Macro9rel <- as.data.frame(make_relative(macro9matrix))
# Create row with maximum relative abundance for each species
colMax <- function(Macro9rel) sapply(Macro9rel, max, na.rm = TRUE)
MaxRel <- colMax(Macro9rel)
Macro9rel <- rbind(Macro9rel, MaxRel = MaxRel)
# Order columns based on the maximum relative abundance
Macro9rel <- Macro9rel[, order(-Macro9rel[nrow(Macro9rel), ])]
# Create 'Other' column with species always less than 2% relative
# abundance
Macro9rel$Others <- rowSums(Macro9rel[, c(17:41)])
MacTax <- Macro9rel[1:9, -c(17:41)]
MacTax$Sample <- rownames(MacTax)
# Create long table
MacTaxLong <- melt(MacTax, id = "Sample")
# Display taxa and samples in the desired order
MacTaxLong$variable <- factor(MacTaxLong$variable, levels = c("Lepetodrilus_fucensis",
"Depressigyra_globulus", "Paralvinella_palmiformis", "Amphisamytha_carldarei",
"Paralvinella_sulfincola", "Provanna_variabilis", "Hydrozoa_non_ID",
"Euphilomedes_climax", "Protomystides_verenae", "Fucaria_striata",
"Munidopsis_alvisca", "Branchinotogluma_tunnicliffeae", "Pardalisca_endeavouri",
"Helicoradomenia_juani", "Clypeosectus_curvus", "Paralvinella_pandorae",
"Others"))
MacTaxLong$Sample <- factor(MacTaxLong$Sample, levels = c("MVw13", "ECw11",
"EMw3", "EMw5", "MVw12", "ECw10", "EMw1", "EMw4", "EMw8"))
MACROcolors17 <- c("#fcd0a1", "#ad2929", "#0fa3b1", "#81f5f5", "#ffea51",
"#d68fd6", "#6144bf", "#f6ae2d", "#8c1c13", "#53917e", "#f26326", "#fffcc9",
"#7098d1", "#cfea9f", "#d26280", "#0686d6", "#979898")
ggplot(MacTaxLong, aes(x = Sample, y = value, fill = variable, order = Sample)) +
geom_bar(stat = "identity", width = 0.9) + scale_fill_manual(values = MACROcolors17,
label = c("L. fucensis", "D. globulus", "P. palmiformis", "A. carldarei",
"P. sulfincola", "P. variabilis", "Hydrozoa non ID", "E. climax",
"P. verenae", "F. striata", "M. alvisca", "B. tunnicliffeae", "P. endeavouri",
"H. juani", "C. curvus", "P. pandorae", "Others")) + theme(axis.text.x = element_text(angle = 90,
hjust = 1, vjust = 0.5), legend.text = element_text(size = 12), legend.title = element_blank()) +
scale_y_continuous(expand = c(0, 0)) + labs(x = NULL, y = "Relative Abundance")
5.2 Meiofauna
5.3 Bacteria
For all microbial plots, the group “Others” includes taxa that are always less than 1% relative abundance. Microbial OTUs were aggregated by high level taxonomies in Excel.
TaxaB23 <- read.csv("Bact23_high_level_taxa.csv", header = TRUE, sep = ",")
TaxaB23LONG <- melt(TaxaB23, id = c("Sample", "Type"))
TaxaB23LONG$variable <- factor(TaxaB23LONG$variable, levels = c("Alphaproteobacteria",
"Gammaproteobacteria", "Deltaproteobacteria", "Proteobacteria", "Bacteroidetes",
"Epsilonbacteraeota", "Acidobacteria", "Actinobacteria", "Aquificae",
"Chloroflexi", "Firmicutes", "Nitrospinae", "Patescibacteria", "Others",
"Bacteria"))
# Worms ordered low to high basal temperature
TaxaB23LONG$Sample <- factor(TaxaB23LONG$Sample, levels = c("bkEM", "bkEC",
"bkMV", "EMf3", "ECf11", "EMf2", "ECf9", "MVf12", "EMf7", "EMf1", "EMw6",
"EMw3", "ECw11", "EMw2", "EMw5", "ECw9", "MVw13", "MVw12", "ECw10",
"EMw7", "EMw1", "EMw4", "EMw8"))
ChromatBact15 <- c("#c60609", "#f3722c", "#f8961e", "#f9c74f", "#9fc681",
"#47682f", "#6bc4a9", "#327e67", "#91a9bd", "#577590", "#b1a0f1", "#7a5ee5",
"#5dbdd5", "#695746", "#928477")
labels <- c(Bkgd = "Bkgd", Diffuse = "Diffuse", Tubeworms = "Tubeworm grabs")
ggplot(TaxaB23LONG, aes(x = Sample, y = value, fill = variable, order = variable)) +
geom_bar(stat = "identity", width = 0.9) + facet_grid(. ~ Type, labeller = labeller(Type = labels),
scales = "free", space = "free") + scale_fill_manual(values = ChromatBact15,
name = "Bacteria", label = c("Alphaproteobacteria", "Gammaproteobacteria",
"Deltaproteobacteria", "Other Proteobacteria", "Bacteroidetes",
"Epsilonbacteraeota", "Acidobacteria", "Actinobacteria", "Aquificae",
"Chloroflexi", "Firmicutes", "Nitrospinae", "Patescibacteria",
"Others", "Unclassified Bacteria")) + theme(axis.text.x = element_text(angle = 90,
hjust = 1, vjust = 0.5), legend.text = element_text(size = 9), legend.title = element_text(size = 12)) +
scale_y_continuous(expand = c(0, 0)) + labs(x = NULL, y = "Relative Adundance")
5.4 Eukarya
TaxaE23 <- read.csv("Euk23_high_level_taxa.csv", header = TRUE, sep = ",")
5.5 Archaea
TaxaA22 <- read.csv("Arch22_high_level_taxa.csv", header = TRUE, sep = ",")
5.6 Balanced microbial community
MG <- read.csv("Combined_micro_major_groups.csv", header = TRUE, sep = ",")
MG$Taxa <- factor(MG$Taxa, levels = c("Bacteria-Gammaproteobacteria", "Bacteria-Alphaproteobacteria",
"Bacteria-Bacteroidetes", "Bacteria-Epsilonbacteraeota", "Bacteria-Deltaproteobacteria",
"Archaea-Thaumarchaeota", "Archaea-Euryarchaeota", "Bacteria-Patescibacteria",
"Bacteria-Chloroflexi", "Eukaryota-Alveolata", "Eukaryota-Rhizaria",
"Bacteria", "Bacteria-Actinobacteria", "Bacteria-Proteobacteria", "Bacteria-Betaproteobacteriales",
"Eukaryota-Opisthokonta", "Eukaryota", "Archaea-Asgardaeota", "Bacteria-Acidobacteria",
"Eukaryota-Amoebozoa", "Bacteria-Aquificae", "Archaea-Crenarchaeota",
"Bacteria-Other", "Archaea-Other", "Eukaryota-Other"))
# Ordered according to hierarchichal clustering dendrogram
MG$source <- factor(MG$source, levels = c("EMf11", "bkMV", "bkEC", "bkEM",
"EMf3", "EMf2", "ECf9", "EMf1", "EMf7", "MVf12", "MVw12", "EMw1", "ECw10",
"EMw7", "EMw4", "EMw8", "MVw13", "ECw11", "ECw9", "EMw3", "EMw6", "EMw2",
"EMw5"))
TAXAcolors25 <- c("#2bd9fe", "#fffcc9", "#45cb85", "#ad2929", "#46237a",
"#f26326", "#0686d6", "#fcd0a1", "#53917e", "#470ff4", "#cc2084", "#63535b",
"#e1dd8f", "#37a4bf", "#cfea9f", "#750d37", "#6144bf", "#7098d1", "#f6ae2d",
"#81f5f5", "#5c5fe8", "#f49390", "#979898", "#a0acad", "#bcb9b8")
ggplot(MG, aes(x = source, y = relabund, fill = Taxa, order = source)) +
geom_bar(stat = "identity", width = 0.9) + scale_fill_manual(values = TAXAcolors25,
labels = c("Gammaproteobacteria", "Alphaproteobacteria", "Bacteroidetes",
"Epsilonbacteraeota", "Deltaproteobacteria", "Thaumarchaeota",
"Euryarchaeota", "Patescibacteria", "Chloroflexi", "Alveolata",
"Rhizaria", "Uncl. Bacteria", "Actinobacteria", "Proteobacteria",
"Betaproteobacteria", "Opisthokonta", "Uncl. Eukaryota", "Asgardaeota",
"Acidobacteria", "Amoebozoa", "Aquificae", "Crenarchaeota", "Other Bacteria",
"Other Archaea", "Other Eukarya")) + theme(axis.text.x = element_text(angle = 90,
hjust = 1, vjust = 0.5), legend.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 9)) + scale_y_continuous(expand = c(0,
0)) + labs(x = NULL, y = "Relative Abundance")
5.7 Details of top 4 bacterial groups
Code only shown for Alphaproteobacteria.
alphas <- read.csv("Alphas.csv", header = TRUE, sep = ",")
alphasLONG <- melt(alphas, id = c("Sample"))
AlphaColors <- c("#cc2084", "#63535b", "#e1dd8f", "#37a4bf", "#cfea9f",
"#750d37", "#6144bf", "#7098d1", "#f6ae2d", "#81f5f5", "#5c5fe8", "#f49390",
"#934b00", "#ffea51", "#488243", "#fffcc9", "#45cb85", "#ad2929", "#46237a",
"#f26326", "#0686d6", "#fcd0a1", "#53917e", "#470ff4", "#CCC9CF", "#CCC9CF",
"#CCC9CF")
alphasLONG$variable <- factor(alphasLONG$variable, levels = c("Rhodobacteraceae",
"SAR11_clade", "unknown_Alphaproteobacteria", "Roseobacter_clade_NAC11_7_lineage",
"Parvibaculales_PS1_clade", "Kordiimonadales", "Robiginitomaculum",
"Euryhalocaulis", "Sedimentitalea", "Bradyrhizobium", "Pseudahrensia",
"Rhodospirillales_AEGEAN_169_marine_group", "Devosiaceae", "Rhizobiales",
"Hellea", "Methylobacterium", "Beijerinckiaceae", "Aestuariispira",
"Nisaeaceae_OM75_clade", "Kiloniellaceae", "Parvularculaceae", "Magnetospiraceae",
"Parvibaculaceae", "Other_Alphaproteobacteria"))
alphasLONG$Sample <- factor(alphasLONG$Sample, levels = c("ECf11", "bkMV",
"bkEC", "bkEM", "EMf3", "EMf2", "ECf9", "EMf1", "EMf7", "MVf12", "MVw12",
"EMw1", "ECw10", "EMw7", "EMw4", "EMw8", "MVw13", "ECw11", "ECw9",
"EMw3", "EMw6", "EMw2", "EMw5"))
ggplot(alphasLONG, aes(x = Sample, y = value, fill = variable, order = variable)) +
geom_bar(stat = "identity", width = 0.9) + scale_fill_manual(values = AlphaColors,
name = NULL, labels = c("Rhodobacteraceae", "SAR11 clade", "Unknown Alphaproteobacteria",
"Roseobacter clade NAC11-7", "Parvibaculales clade PS1", "Kordiimonadales",
"Robiginitomaculum", "Euryhalocaulis", "Sedimentitalea", "Bradyrhizobium",
"Pseudahrensia", "Rhodospirillales clade AEGEAN-169", "Devosiaceae",
"Rhizobiales", "Hellea", "Methylobacterium", "Beijerinckiaceae",
"Aestuariispira", "Nisaeaceae clade OM75", "Kiloniellaceae", "Parvularculaceae",
"Magnetospiraceae", "Parvibaculaceae", "Other Alphaproteobacteria")) +
theme(axis.text.x = element_text(size = 11, angle = 90, hjust = 1,
vjust = 0.5), legend.text = element_text(size = 9)) + labs(x = NULL,
y = "Relative Abundance")
gammas <- read.csv("Gammas.csv", header = TRUE, sep = ",")
bacteroidetes <- read.csv("Bacteroidetes.csv", header = TRUE, sep = ",")
epsilons <- read.csv("Epsilons.csv", header = TRUE, sep = ",")
6 Enriched species/OTUs
Species and OTUs that differentiated between different sample types (tubeworm grabs, diffuse fluids and background fluids) and grab sample temperature categories (highT, lowT) were identified using ALDEx2. Microbes were analyzed from all 13 grab samples (6 highT, 7 lowT) and 10 fluid samples (7 diffuse, 3 background) using both individual domains and the QPCR-balanced microbial assemblage. Fauna were analyzed from 9 grab samples (5 highT, 4 lowT). We include an analysis example below, but do not include all tests for brevity.
library(ALDEx2)
# create the 2-state condition for comparison
tempx2 <- as.vector(metadata9$T_category)
tempx2## [1] "low" "low" "high" "high" "high" "high" "low" "high" "low"# ALDEx2 modular mode
flip_mac9 <- as.data.frame(t(macro9))
clr_mac9byT <- aldex.clr(flip_mac9, tempx2, mc.samples = 128, denom = "all")
kw_mac9byT <- aldex.kw(clr_mac9byT)
kw_mac9byT.effect <- aldex.effect(clr_mac9byT, tempx2, verbose = FALSE)
mac9byT.all <- data.frame(kw_mac9byT, kw_mac9byT.effect)
# mac9byT.all has details of differential taxa between highT and lowT
# (plotted in red below)
aldex.plot(mac9byT.all, type = "MW", test = "glm", xlab = "Dispersion",
ylab = "Difference", all.cex = 1.5, all.pch = 16, all.col = "lightgrey",
rare.col = "black")
- In this case, macrofauna enriched in highT assemblages (effect < -1, Benjamini-Hochberg corrected p-values < 0.1 for GLM test):
- Paralvinella sulfincola
- Paralvinella pandorae
- Depressigyra globulus
- Paralvinella palmiformis
- Branchinotogluma tunnicliffeae
- Macrofauna enriched in lowT assemblages (effect > 1):
- Amphisamytha caldarei
- Helicoradomenia juani
7 Plotting enriched species/OTUs
Microbes presented here are from ALDEx2 tests using the QPCR-balanced microbial assemblage. Taxa with absolute effect values >1 were compiled in Excel for plotting.
7.1 HighT vs. lowT grabs
htlteBAL <- read.csv("enriched_HTvsLT_grabs.csv", header = TRUE, sep = ",")
htlteBAL$range <- ifelse(htlteBAL$Effect < 0, "below", "above") # 'above'= low T 'below'= high T
htlteBAL <- htlteBAL[order(htlteBAL$Effect), ]
htlteBAL$Identity <- factor(htlteBAL$Identity, levels = htlteBAL$Identity) # keeps the order for the plot rather than going alphabetical
ggplot(htlteBAL, aes(x = Identity, y = Effect, label = Effect)) + geom_bar(stat = "identity",
aes(fill = Domain), width = 0.6) + scale_fill_manual(name = NULL, labels = c("Archaea",
"Bacteria", "Eukarya", "Macrofauna", "Meiofauna"), values = c(Archaea = "#abdafc",
Bacteria = "#0686d6", Eukarya = "#cfea9f", Macrofauna = "#f6ae2d",
Meiofauna = "#f26326")) + theme(axis.text.x = element_text(size = 8),
axis.title.x = element_text(size = 9), axis.text.y = element_text(size = 6),
axis.title.y = element_blank(), legend.position = "none") + ggtitle("HighT vs. LowT grabs") +
scale_x_discrete(position = "bottom") + coord_flip()
7.2 HighT grabs vs. diffuse
All 6 highT grab samples and 3 associated diffuse fluid samples.
7.3 LowT grabs vs. diffuse
All 7 lowT grab samples and 4 associated diffuse fluid samples.
8 Congruence
Symmetric procrustes analysis (function ‘procrustes’) and significance testing (function ‘protest’) performed on pairs of NMDS ordinations.
8.1 Eight sample NMDS
Archaeal reads were insufficient from sample ECw11 so NMDS plots without this sample for the other size classes and microbial domains will be necessary for comparisons.
8.2 Nine sample, QPCR-balanced microbial NMDS
Envfit.
##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## Temperature 0.91082 0.41280 0.4431 0.026 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## NMDS1 NMDS2
## SourceTubeworms 0.0000 0.0000
## Sample_typebk -25.9926 -0.0094
## Sample_typed 51.9852 0.0189
## Vent_fieldClam_Bed 51.9647 0.0370
## Vent_fieldMain_Endeavour -41.5885 0.0105
## Vent_fieldMiddle_Valley 52.0066 -0.0633
## T_categoryhigh 51.9988 -0.0145
## T_categorylow -64.9985 0.0181
## Substratumbasalt -65.0383 -67.4684
## Substratumsulfide 18.5824 19.2767
##
## Goodness of fit:
## r2 Pr(>r)
## Source 0.0000 1.000
## Sample_type 0.1000 0.469
## Vent_field 0.1600 0.845
## T_category 0.2502 0.025 *
## Substratum 0.1857 0.190
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 9998.3 Macrofauna vs. Meiofauna
library(vegan)
clr_mac9meio9NMDS <- procrustes(clr_mac9.NMDS, clr_meio9.NMDS, scale = FALSE,
choices = c(1, 2), symmetric = TRUE)
par(mfrow = c(1, 2))
plot(clr_mac9meio9NMDS, kind = 1)
plot(clr_mac9meio9NMDS, kind = 2)
Testing for significance.
protest(clr_mac9.NMDS, clr_meio9.NMDS, scores = "sites", permutations = how(nperm = 999))##
## Call:
## protest(X = clr_mac9.NMDS, Y = clr_meio9.NMDS, scores = "sites", permutations = how(nperm = 999))
##
## Procrustes Sum of Squares (m12 squared): 0.2588
## Correlation in a symmetric Procrustes rotation: 0.8609
## Significance: 0.005
##
## Permutation: free
## Number of permutations: 9998.4 Macrofauna vs. Bacteria
##
## Call:
## protest(X = clr_mac9.NMDS, Y = clr_bact9.NMDS, scores = "sites", permutations = how(nperm = 999))
##
## Procrustes Sum of Squares (m12 squared): 0.6513
## Correlation in a symmetric Procrustes rotation: 0.5905
## Significance: 0.096
##
## Permutation: free
## Number of permutations: 9998.5 Macrofauna vs. Archaea
##
## Call:
## protest(X = clr_mac8.NMDS, Y = clr_arch8.NMDS, scores = "sites", permutations = how(nperm = 999))
##
## Procrustes Sum of Squares (m12 squared): 0.6865
## Correlation in a symmetric Procrustes rotation: 0.5599
## Significance: 0.02
##
## Permutation: free
## Number of permutations: 9998.6 Macrofauna vs. Eukarya
##
## Call:
## protest(X = clr_mac9.NMDS, Y = clr_euk9.NMDS, scores = "sites", permutations = how(nperm = 999))
##
## Procrustes Sum of Squares (m12 squared): 0.9062
## Correlation in a symmetric Procrustes rotation: 0.3062
## Significance: 0.818
##
## Permutation: free
## Number of permutations: 9998.7 Meiofauna vs. Bacteria
##
## Call:
## protest(X = clr_meio9.NMDS, Y = clr_bact9.NMDS, scores = "sites", permutations = how(nperm = 999))
##
## Procrustes Sum of Squares (m12 squared): 0.6047
## Correlation in a symmetric Procrustes rotation: 0.6287
## Significance: 0.057
##
## Permutation: free
## Number of permutations: 9998.8 Meiofauna vs. Archaea
##
## Call:
## protest(X = clr_meio8.NMDS, Y = clr_arch8.NMDS, scores = "sites", permutations = how(nperm = 999))
##
## Procrustes Sum of Squares (m12 squared): 0.8721
## Correlation in a symmetric Procrustes rotation: 0.3576
## Significance: 0.285
##
## Permutation: free
## Number of permutations: 9998.9 Meiofauma vs. Eukarya
##
## Call:
## protest(X = clr_meio9.NMDS, Y = clr_euk9.NMDS, scores = "sites", permutations = how(nperm = 999))
##
## Procrustes Sum of Squares (m12 squared): 0.776
## Correlation in a symmetric Procrustes rotation: 0.4733
## Significance: 0.268
##
## Permutation: free
## Number of permutations: 9998.10 Bacteria vs. Eukarya
##
## Call:
## protest(X = clr_bact9.NMDS, Y = clr_euk9.NMDS, scores = "sites", permutations = how(nperm = 999))
##
## Procrustes Sum of Squares (m12 squared): 0.3398
## Correlation in a symmetric Procrustes rotation: 0.8125
## Significance: 0.001
##
## Permutation: free
## Number of permutations: 9998.11 Archaea vs. Eukarya
##
## Call:
## protest(X = clr_arch8.NMDS, Y = clr_euk8.NMDS, scores = "sites", permutations = how(nperm = 999))
##
## Procrustes Sum of Squares (m12 squared): 0.9796
## Correlation in a symmetric Procrustes rotation: 0.1429
## Significance: 0.378
##
## Permutation: free
## Number of permutations: 9998.12 Bacteria vs. Archaea
##
## Call:
## protest(X = clr_bact8.NMDS, Y = clr_arch8.NMDS, scores = "sites", permutations = how(nperm = 999))
##
## Procrustes Sum of Squares (m12 squared): 0.9622
## Correlation in a symmetric Procrustes rotation: 0.1944
## Significance: 0.756
##
## Permutation: free
## Number of permutations: 9998.13 Microbes vs. Macrofauna
##
## Call:
## protest(X = clr_mac9.NMDS, Y = clr_balMic9.NMDS, scores = "sites", permutations = how(nperm = 999))
##
## Procrustes Sum of Squares (m12 squared): 0.8275
## Correlation in a symmetric Procrustes rotation: 0.4153
## Significance: 0.445
##
## Permutation: free
## Number of permutations: 9998.14 Microbes vs. Meiofauna
##
## Call:
## protest(X = clr_meio9.NMDS, Y = clr_balMic9.NMDS, scores = "sites", permutations = how(nperm = 999))
##
## Procrustes Sum of Squares (m12 squared): 0.7268
## Correlation in a symmetric Procrustes rotation: 0.5226
## Significance: 0.093
##
## Permutation: free
## Number of permutations: 9999 Proportionality
Microbial OTUs were binned by shared taxonomic assignment into 719 unique taxonomic identities. Prior to proportionality, rare taxa (i.e. those with counts of <2 in 7+ samples) were removed to minimize false discovery rate, resulting in final numbers of taxa for analysis of 335 unique microbial taxa, 14 macrofauna, and 17 meiofauna. Only rho values >0.75 are used in Cytoscape network mapping.
library(propr)
NoRares <- read.csv("Non_rares_for_propr.csv", header = T, row.names = 1)
# Imputation of zeros
NoRares_cmultRepl <- cmultRepl(NoRares, method = "CZM", output = "p-counts")## No. corrected values: 85# Centered log-ratio transformation
clr_NoRares <- apply(NoRares_cmultRepl, 2, function(x) {
log(x) - mean(log(x))
})
# Proportionality analysis
NoRares.propr <- propr(counts = NoRares_cmultRepl, metric = "rho", p = 1000)
NoRaresCorrALL <- getResults(NoRares.propr)10 Session Info
sessionInfo()## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ALDEx2_1.20.0 propr_4.2.6 zCompositions_1.3.4
## [4] truncnorm_1.0-8 NADA_1.6-1.1 survival_3.1-12
## [7] MASS_7.3-51.6 compositions_2.0-0 funrar_1.4.1
## [10] ggpubr_0.4.0 ggrepel_0.8.2 ggplot2_3.3.2
## [13] reshape2_1.4.4 vegan_2.5-6 lattice_0.20-41
## [16] permute_0.9-5 knitr_1.29
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-148 bitops_1.0-6
## [3] matrixStats_0.56.0 GenomeInfoDb_1.24.2
## [5] tensorA_0.36.1 tools_4.0.2
## [7] backports_1.1.9 R6_2.4.1
## [9] BiocGenerics_0.34.0 mgcv_1.8-31
## [11] colorspace_1.4-1 withr_2.2.0
## [13] gridExtra_2.3 tidyselect_1.1.0
## [15] bayesm_3.1-4 curl_4.3
## [17] compiler_4.0.2 Biobase_2.48.0
## [19] formatR_1.7 DelayedArray_0.14.1
## [21] labeling_0.3 scales_1.1.1
## [23] DEoptimR_1.0-8 robustbase_0.93-6
## [25] stringr_1.4.0 digest_0.6.25
## [27] foreign_0.8-80 rmarkdown_2.3
## [29] rio_0.5.16 XVector_0.28.0
## [31] pkgconfig_2.0.3 htmltools_0.5.0
## [33] rlang_0.4.7 readxl_1.3.1
## [35] generics_0.0.2 farver_2.0.3
## [37] BiocParallel_1.22.0 dplyr_1.0.2
## [39] zip_2.1.1 car_3.0-9
## [41] RCurl_1.98-1.2 magrittr_1.5
## [43] GenomeInfoDbData_1.2.3 Matrix_1.2-18
## [45] Rcpp_1.0.5 munsell_0.5.0
## [47] S4Vectors_0.26.1 abind_1.4-5
## [49] lifecycle_0.2.0 stringi_1.5.3
## [51] yaml_2.2.1 carData_3.0-4
## [53] SummarizedExperiment_1.18.2 zlibbioc_1.34.0
## [55] plyr_1.8.6 grid_4.0.2
## [57] parallel_4.0.2 forcats_0.5.0
## [59] crayon_1.3.4 haven_2.3.1
## [61] cowplot_1.1.0 splines_4.0.2
## [63] hms_0.5.3 pillar_1.4.6
## [65] GenomicRanges_1.40.0 ggsignif_0.6.0
## [67] stats4_4.0.2 glue_1.4.2
## [69] evaluate_0.14 data.table_1.13.0
## [71] vctrs_0.3.4 cellranger_1.1.0
## [73] gtable_0.3.0 purrr_0.3.4
## [75] tidyr_1.1.2 xfun_0.17
## [77] openxlsx_4.1.5 broom_0.7.0
## [79] rstatix_0.6.0 tibble_3.0.3
## [81] IRanges_2.22.2 cluster_2.1.0
## [83] ellipsis_0.3.1